package cmlabs.numbertheoryutils.algorithms;

import java.math.BigInteger;
import java.util.ArrayList;
import java.util.List;
import java.util.Random;
import java.util.TreeMap;

/**
 * Created by IntelliJ IDEA.
 * User: jte
 * Date: 9.4.11
 * Time: 18.24
 * To change this template use File | Settings | File Templates.
 */
public class Factorization {

    public static TreeMap<BigInteger, Integer> factorize(BigInteger n) {
        factors.clear();
        makeFactorization(n);
        return (TreeMap<BigInteger, Integer>) factors.clone();
    }

    private static void makeFactorization(BigInteger n) {
        if (!n.equals(BigInteger.ONE)) {
            if (n.isProbablePrime(20)) {
                if (factors.containsKey(n)) {
                    factors.put(n, factors.get(n) + 1);
                } else {
                    factors.put(n, 1);
                }
            } else {
                BigInteger divisor = n.and(BigInteger.ONE);
                if (divisor.equals(BigInteger.ZERO))
                    divisor = BigInteger.valueOf(2);
                if (divisor.equals(BigInteger.ONE))
                    divisor = pollardMonteCarlo(n);
                if (divisor.equals(BigInteger.ONE))
                    divisor = pollardRHO(n);
                if (divisor.equals(BigInteger.ONE))
                    divisor = pollardBent(n);
                makeFactorization(divisor);
                makeFactorization(n.divide(divisor));
            }
        }
    }

    private static BigInteger pollardMonteCarlo(BigInteger n) {
        BigInteger m = BigInteger.valueOf(100);
        BigInteger b = randomBigInteger(n.bitLength()).mod(m.subtract(BigInteger.valueOf(2))).add(BigInteger.valueOf(2));
        List<BigInteger> primes = new ArrayList<BigInteger>();
        BigInteger mMax = m.multiply(BigInteger.ONE);
        primes.add(BigInteger.valueOf(3));
        for (BigInteger prime = BigInteger.valueOf(5); prime.compareTo(m) <= 0; prime = prime.add(BigInteger.ONE)) {
            if (prime.isProbablePrime(20)) {
                primes.add(prime);
            }
        }
        BigInteger g = BigInteger.ONE;
        for (int i = 0; i < primes.size() && g.equals(BigInteger.ONE); ++i) {
            BigInteger current = primes.get(i);
            while (current.compareTo(n) <= 0) {
                current = current.multiply(primes.get(i));
            }
            current = current.divide(primes.get(i));
            b = b.modPow(current, n);
            g = b.subtract(BigInteger.ONE).abs().gcd(n);
            if (g.equals(n)) {
                g = BigInteger.ONE;
            }
        }
        return g;
    }

    private static BigInteger pollardBent(BigInteger n) {
        int iterations = 33;
        BigInteger b0 = randomBigInteger(n.bitLength()).mod(n), b1 = b0.multiply(b0).add(BigInteger.valueOf(2)).mod(n), a = b1.multiply(BigInteger.ONE);
        for (int i = 0, seriesLength = 1; i < iterations; ++i, seriesLength <<= 1) {
            BigInteger g = b1.subtract(b0).abs().gcd(n);
            for (int j = 0; j < seriesLength && !g.equals(n) && !g.equals(BigInteger.ONE); ++j) {
                b1 = b1.multiply(b1).add(BigInteger.valueOf(2)).mod(n);
                g = b1.subtract(b0).abs().gcd(n);
            }
            b0 = a.multiply(BigInteger.ONE);
            a = b1.multiply(BigInteger.ONE);
            if (!g.equals(BigInteger.ONE) && !g.equals(n)) {
                return g;
            }
        }
        return BigInteger.ONE;
    }

    private static BigInteger pollardRHO(BigInteger n) {
        int iterations = 10000;
        BigInteger b0 = randomBigInteger(n.bitLength()).mod(n), b1 = b0.add(BigInteger.ZERO), g;
        b1 = b1.pow(2).mod(n).add(BigInteger.ONE).mod(n);
        g = b1.subtract(b0).abs().gcd(n);
        for (int i = 0; i < iterations; ++i) {
            b0 = b0.pow(2).mod(n).add(BigInteger.ONE).mod(n);
            b1 = b1.pow(2).mod(n).add(BigInteger.ONE).mod(n);
            g = b1.subtract(b0).abs().gcd(n);
        }
        return g;
    }

    private static BigInteger randomBigInteger(int n) {
        byte[] temporary = new byte[n / 8 + 1];
        new Random(System.currentTimeMillis()).nextBytes(temporary);
        return new BigInteger(temporary);
    }

    private static TreeMap<BigInteger, Integer> factors = new TreeMap<BigInteger, Integer>();
}
